function [pdf,logpdf]=invwishlub(sigma,d,S) 
%function [pdf,logpdf]=invwishlub(sigma,d,S) 
% INVWISHLUB = Inverse Wishart p.d.f.
% sigma = positive definite
% S = q x q symmetric
% ( Corresponds to d in Lubrano et. al)
% pdf = probability density function, properly normalized
% Formula Taken from Lubrano et. al page 305 
% NOTE: For the unidimensional case this transaltes to 
% (1/sigma) ~ IG ( d/2,2/S)
% SEE also invwpdf.m 
% Double checked 6/18/03 
q=size(sigma,1); 
vet=((d+1 -[1:q]))/2;
logk=zeros(4,1); 
logk(1)=d*q*log(2)/2; 
logk(2)=q*(q-1)*log(pi)/4; 
logk(3)=sum(gammaln(vet)); 
logk(4)=(-d/2)*log(det(S)); 
logk=sum(logk);
logexpo=-trace(S/sigma)/2; 
logdeter=-(d+q+1)*log(det(sigma))/2; 
logpdf = logexpo +logdeter - logk; 
pdf=exp(logpdf); 
